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ABSTRACT 

We describe a simple method for estimating the vertical column density i n Smoothed 
Particle Hy drodynamics (SPH) simulations of discs. As in the method of Stamatellos 



1 INTRODUCTION 



Smooth Particle Hydrodynamics (SPH) (Lucy 1977 Gin- 
gold k, Monaghan|1977 ) is a Lagrangian technique for simu- 
lating fluid flows using a particle representation. This tech- 
nique assigns each particle a mass, position and internal en- 
ergy and then interpolates state variables, such as density, by 
"smoothing" over neighbouring particles. Gravity has been 
incorporated into the SPH formalism, allowing SPH to be 
used to simulate astrophysical fluids. However, the thermal 
evolution of many astrophysical systems is often governed 
by radiative transfer effects, in addition to hydrodynamic 
energy transfer. As an accurate description of a system's 
thermal evolution is of vital importance in many astrophys- 
ical systems (e.g. accretion discs, stellar gas clouds) a lot of 
effort has recently been made to add radiative transfer to 



et al. (2007), the column density is estimated using pre-computed local quantities and 
is then used to estimate the radiative cooling rate. The cooling rate is a quantity of 
considerable importance, for example, in assessing the probability of disc fragmenta- 
tion. Our method has three steps: (i) the column density from the particle to the mid 
plane is estimated using the vertical component of the gravitational acceleration, (ii) 
the "total surface density" from the mid plane to the surface of the disc is calculated, 
(iii) the column density from each particle to the surface is calculated from the dif- 
ference between (i) and (ii). This method is shown to greatly improve the accuracy of 
column density estimates in disc geometry compared with the method of Stamatellos. 
On the other hand, although the accuracy of our method is still acceptable in the case 
of high density fragments formed within discs, we find that the Stamatellos method 
performs better than our method in this regime. Thus, a hybrid method (where the 
method is switched in regions of large over-density) may be optimal. 

underlying physics, it has a very low computational cost 
and so simulations including approximating radiative cool- 
ing can be run without sacrificing spatial or temporal reso- 
lution. Recently, Stamatellos et al have proposed a method 
that improves upon the cooling time prescription without 
significantly increasing the computational cost (Stamatellos] 
2007| . Forgan et al extended this method to include 



et al. 



SPH (Oxley Woolfson 2003 Whitehouse Bate 2004a 



Stamatellos et al.||2007| [Forgan et al.||2009| ). Unfortunately, 
a full three dimensional, frequency dependent description of 
radiative transfer is currently computationally impossible, 
so simplifying assumptions have to be made. 

Since SPH is a Lagrangian method, the purpose of mod- 
eling radiative transfer effects is to provide a net cooling (or 
heating) rate per particle as a result of radiative processes. 
For example, a popular choice for simulating self- gravitating 
discs is to use a cooling time prescription, where the radia- 
tive cooling rate U is set equal to where U is the cooling 

''cooZ 

rate per unit mass and tcooi is simply parameterised as a pre- 
scribed multiple of the local dynamical timescale (Rice et al. 



heat transfer between particles, by combining the Stamatel- 
los method with the flux limited diffusion (FLD) method 
dForgan et al.||2009j ). 

The Stamatellos method estimates the optical depth by 
assuming that the relationship between the two locally com- 
puted variables, density and gravitational potential, and the 
optical depth is the same as it is for a mass- weighted average 
of that relationship over a self- gravitating poly tropic sphere. 
They then estimate the local cooling rate using only the local 
temperature and optical depth. This estimate tends to the 
radiative diffusion approximation at high optical depths but 
does not involve calculating noisy derivatives as is required 



in a proper implementation of radiative diffusion (White- 
house & Bate 2004b). The aim of such methods is not to 



2003 ) . Although such a description grossly oversimplifies the 



achieve perfect agreement with the full radiative transfer 
equations, but to achieve an accuracy that is comparable 
in magnitude to the other uncertainties such as those as- 
sociated with the grey approximation and the appropriate 
values of the frequency averaged opacity. 

However, although the Stamatellos method has proved 
successful for modelling spherical systems, Wilkins & Clarke 
have shown that it can systematically underestimate the 
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cooling in disc geometries by as much as a factor of four 
in the region of the mid plane, where column density es- 



timation is critical to estimating the cooling (jWilkins & 
Clarke|2012 ). This underestimate stems from the Stamatel- 
los method overestimating the column density, S, which ap- 
pears as a quadratic term in the expression for the cooling 
rate in the optically thick limit. An example of where radia- 
tive transfer is important to the evolution of a disc system 
is the study of gravitational instabilities in proto-planetary 
disc. The cooling time prescription has been used to study 
these discs in great detail, but a more accurate description of 
the cooling in such systems is necessary to improve our un- 



derstanding of these gravitational instabilities ( Forgan et al. 



2011). 



In this paper we propose a variation on the Stamatellos 
method for calculating the cooling rate for the special case 
of disc geometries, by improving the estimate of the column 
density in this case. Our method still only requires the use of 
local quantities, already calculated by the gravitational and 
hydrodynamic codes and as such remains computationally 
inexpensive. The paper is organized as follows. In section |2] 
we describe our new method in detail. In section [3] we test 
our method on a series of discs for which all quantities of 
interest can be obtained analytically or semi-analytically. In 
section |4] we test our method on realistic disc simulations 
that have been evolved long enough to either fragment or 
reach marginal stability. Finally, section |5] summarises our 
conclusions. 



2 THE METHOD 

For a typical astrophysical disc, most of the cooling occurs 
at the disc's surface. Therefore, we assume that the greatest 
contribution to the radiative cooling of a particle within such 
a disc comes from energy radiated vertically out of the disc, 
rather than along the optically thick mid-plane. Assuming 
this is true, we estimate the cooling rate of each particle 
using the optical depth, r, along a vertical path from the 
particle to the disc's surface. The true optical depth is given 
by T = J Hi{z)p{z)dz. In order to avoid evaluating this costly 
integral we make the approximation that t — Hi j pdz — kT^, 
where S is the column density to the surface and k is the 
local opacity. Next, we follow Stamatellos et al in defining 
the cooling rate per unit mass to be: 



(1) 



where, r is the optical depth to the surface of the disc, a 
is the Stefan Boltzman constant, S is the column density 
from the particle to the surface of the disc and /? is the 
local Planck-mean opacity. The To term is a background 
temperature below which particles are not allowed to cool. 

In the regime where > we find that there are 
two limiting cases. The first limiting case is the optically 
thin limit, where rS ^ At"^. In this limit, the cooling simply 
reduces to 



(2) 



which is just the cooling rate for an isolated particle, in an 
environment with opacity n. However, note that this is not a 
good approximation to the cooling rate per unit mass for an 



optically thin layer on top of an optically thick disc, which is 
a commonly encountered situation. Fortunately, in practice 
the thermodynamics of such layers are often controlled by 
the background temperature, Tq. 

In the optically thick limit, rS ^ and equation [l] 
becomes 



U ^ 



St 



(3) 



which is an approximation to the commonly used diffusion 



approximation (see Mihalas ( 1970|, section 2.3). To see why, 
consider the optically thick limit, in which the radiative flux 
is given by: 



F = 



and the corresponding cooling rate per unit mass as: 



-V.F 



(4) 



(5) 



If integration and differentiation over z are replaced by 
simple multiplication and division by an effective vertical 
scale height, H, then an expression of the same form as [3] is 
obtained. 



U = -V.F! 
P 



(6) 



Note that because U oc any inaccuracy in the col- 
umn density to the surface, S, will have a large effect on the 
inaccuracy of the cooling rate. 

S is calculated in a three step process. Firstly, an es- 
timate of the column density between the particle and the 
mid-plane is obtained. This is done using only the verti- 
cal component of the gravitational acceleration, which is 
already calculated by the simulation making this step es- 
sentially "free" from a computational standpoint. Secondly, 
a total surface density map (i.e. map of column density from 
the mid-plane to the surface) is computed. Finally, the col- 
umn density to the surface is calculated by subtracting the 
column density to the mid-plane from the total surface den- 
sity. 



2.1 Estimating the column density to the 
mid-plane 

We estimate the column density of each particle to the mid- 
plane using the approximation: 



9z 



GMz 



■ AttGTi' sign{z) 



(7) 



where M is the mass of the central star, r — x'^ -\- y'^ -\- z'^ 
and is the column density between the point at height 
z and the mid-plane. The first term in this expression is 
the contribution from the central star, whereas the second 
term comes from assuming that the disc is infinite in extent 
and has the same vertical density structure everywhere. In 
reality, there will be further corrections to this expression 
due to radial variation in the disc's vertical density structure 



and its finite extent (see e.g. Appendix of Bertin & Lodato 
( 1999j ).) However, we expect these contributions to be small 
in most cases, an assumption that will be investigated in 
detail in section [S] 
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As the gravitational acceleration, Qz^ is already calcu- 
lated by the simulation for every particle, we can re-arrange 
equation ([7| to estimate the column density from the parti- 
cle to the mid-plan^ 



E' 



sign{z 
4nG 



( ^ ^^^^ 



(8) 



which must be performed by any SPH code. As such, even in 
the worst case scenario this step will be no more expensive 
than the 3D density calculation already performed by the 
code. Furthermore, as there is little gain in calculating the 
surface density to greater accuracy than the other sources 
of error in our method (see discussion below), a computa- 
tionally cheaper density estimator can usually be used. 



2.2 Constructing a surface density map 

In order to calculate the cooling, we need to know the col- 
umn density between the particle and the surface of the 
disc, and therefore the column density to the mid-plane, es- 
timated in equation ([8| above, has to be subtracted from an 
estimate of the total column density. This map can be cal- 
culated in a number of different ways, discussed at length in 
Appendix [B] and the section below. In this paper, we calcu- 
late the map by projecting all particles onto the mid-plane 
and then interpolating for each particle from a density map 
evaluated at the location of a subset of particles (10%). The 
surface density for the remaining 90% of points are inter- 
polated when needed. Once obtained, we can estimate the 
column density to the surface, S, by: 



(9) 



2.3 Computational efficiency 



The construction of a total surface density map, has the 
potential to be computationally expensive. The potentially 
expensive aspect of this step is the estimation of the sur- 
face density at each particle from the particle positions and 
masses, projected onto the mid plane. There already ex- 
ists an extensive literature devoted to solving this problem, 
which we will not attempt to reproduce in detail here. Im- 
portantly, the choice of method for estimating density in- 
volves a trade off between accuracy, computational speed 



and the ability to recover sharp density gradients (see ( Fer- 
dosi et aT]|2011 ) for a discussion of these trade-offs in the 
context of astronomical datasets). 

SPH density estimation methods give excellent accu- 
racy, but a high computational cost (0(A^^), e.g. the "DED- 
ICA" method in [Ferdosi et al.| ( |2011| ))). At the other ex- 
treme, grid based methods offer excellent computational ef- 
ficiency (0(A^)), but with reduced accuracy (e.g. the "MBE" 
method 



Ferdosi et al. 



( 20TT| )) 



The method we use 



"kNN" method described in Ferdosi et al. (2011), which has 
0{NlogN) scaling. 

The computational cost of any of the above method 
can be reduced, at the cost of "smoothing out" the density 
distribution, by evaluating the density at a subset of points 
and then interpolating the values at the particle locations. 

It is important to note that the calculation of a to- 
tal surface density map is a two-dimensional version of the 
calculation of the three-dimensional density of each particle. 



^ Note that this assumes that the value of gz calculated by the 
code includes the contribution from the star. If it does not, the 
term due to the star's potential should be dropped from equations 
|7|and|8] 



2.4 Sources of inaccuracy 

The accuracy of this method can be affected by a number 
of different factors. Each of these potential sources of error 
will be tested independently in Section [S] 

i) The gravitational force calculation in self-gravitating 
SPH codes does not simply consist of summing pairwise in- 
teractions. In particular, the usual gravitational force is typ- 
ically "softened" for small separations and distant particles 
are grouped together in a 'tree' structure, e.g. an oct-tree 



(Barnes & Hut 1986 ). How aggressively particles are grouped 



and the extent of the gravitational softening is controlled 
by the user. As our method relies on the gravity estimated 
from the tree code, inaccuracies in the gravitational force 
may propagate in our method also. 

ii) The discrete realisation of a given density distribu- 
tion may produce accelerations that differ from those pro- 
duced by a continuous density field since quantising the den- 
sity distribution inevitably introduces "dumpiness", which 
affects the resulting gravitational acceleration. Note that 
this error decreases with increasing resolution. The gravita- 
tional softening mentioned above also mitigates this effect. 

iii) No physical disc is really an infinite slab of con- 
stant height, with density structure depending only on z. 
Deviations from this structure will obviously influence the 
accuracy of our method. 

iv) Any inaccuracy in our map of surface densities Smap 
will lead to inaccuracies in our final estimate of E. This 
inaccuracy will arise from a combination of interpolation 
error (as the map is only calculated for 10% of particles) 
and error due to the calculation of the surface density itself. 

To quantify the size of these relative errors, we note 
that effects i)-iii) produce errors in the estimate of S^, while 
iv) produces errors m Zjrnap- We represent the size of these 
errors as ST>' and Sl^map, respectively. We then see that the 
relative error in our estimate of S is: 



here (describe in appendix |B| is intermediate in complex- § f ^^^^^ \ 

ity between these two extremes and is very similar to the \ E / 



(10) 



Since E is always smallest at large z, it follows that the 
fractional inaccuracy of our estimate is generally worst far 
from the mid-plane. 



3 ANALYTIC TESTS 

As outlined in the previous section, there are a number of 
possible sources of inaccuracy in our column density esti- 
mate. These are: inaccuracies in the gravitational tree code, 
the quantization of a continuous system, the "infinite slab" 
approximation and inaccuracy of the surface density map. 

In this section we construct a series of discs for which 
the gravitational acceleration, density and column density 
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can be calculated analytically or semi-analytically. As de- 
scribed in Appendix [a] these discs have power law density 
profiles between adjustable inner and outer radii, with a 
Gaussian variation of density with z at each radius. The 
scale height of this Gaussian is adjusted so that it corre- 
sponds to the hydrostatic equilibrium profile of a locally 
isothermal non-self gravitating disc where the temperature 
is a power law function of radius. Because we are given a 
continuous, analytic density distribution, we can integrate 
it semi-analytically to calculate the resulting vertical com- 



ponent of the gravitational acceleration (equation A7 ) . Fur- 
thermore, we can analytically determine the gravitational 
acceleration for the continuous density distribution under 
the infinite slab approximation (equation[7|without the term 
due to the star, where is given by equation |A10 ). We also 
calculate the vertical component of gravitational accelera- 
tion using the oct-tree on the discretized distribution. We 
will henceforth refer to the gravitational acceleration calcu- 
lated by integration of the continuous density distribution 
as the "continuous gravitational acceleration" and the value 
calculated using the oct-tree as the "oct-tree gravitational 
acceleration" . Combining these three estimates of Qz , we can 
independently test the sources of error identified in section 
|2.4| We also test the accuracy of our total surface density 
map by comparing it to the analytic value given by |A10| 

We emphasise that these discs do not represent realistic 
physical conditions that would result from the steady state 
evolution of accretion discs, but are designed to test our 
method for calculating the column density. As such, many 
of these discs will have non-uniform values of the Toomre Q 
parameter ( Toomre|1964| ), significantly different from unity 
and would change significantly if allowed to evolve with time. 

The properties of our analytic discs are summarised in 
Table [l] In each case the density structure (Equation Al ) 
corresponds to a hydrostatic structure of uniform temper- 
ature with mid-plane density scaling as R~'^'^. In running 
the simulations, typical parameters for the tree gravity were 
used (adaptive gravitational softening on the SPH smooth- 
ing length /i, Barnes & Hut opening angle threshold of 0.3). 



3.1 Varying scale heights 

We first test how the accuracy of our method depends on 
disc thickness, where thickness is measured by H/R and H 
is the scale height of the disc. We consider a thin disc with 
0.025 < H/R < 0.11 and a thick disc with 0.76 < H/R < 
3.4. 

Figures [l] and [3] show the inaccuracy in the vertical 
component of gravitational acceleration due to making the 



continuous density distribution (equation Al ) discrete and 
demonstrate that, as expected, the errors are highest in thin- 
ner discs, particularly close to the disc mid-plane where any 
given particle is most affected by the "lumpy" nature of its 
environment (see Figures [5] and |4|. The effect is exacerbated 
by the fact that ^ in the densest parts of the disc, so 
the relative errors are largest here. Note that such quanti- 
zation errors are reduced by increasing the resolution (see 
section 



3.2) 



Next we investigate the inaccuracy in the gravitational 
acceleration due to use of the infinite slab approximation. 
Figure [5] show that this error is small for the thin disc, with 
most points being accurate to within 50% of the expected 
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Figure 1. The log2 ratio of the vertical component of the oct- 
tree and continuous gravitational accelerations, as a function of 
cylindrical radius for a subset of 5000 particles in the 'thin disc' 
calculation {H/R = ^m^^R/Ri). The boxplot on the left of the 
plot shows the distribution of errors. The box bounds the central 
50% of particles (the edges are the 25% and 75% percentiles), the 
red line is the median and the whiskers mark the 5% and 95% 
percentiles. 



Spatial distribution of errors due 
to quantization (Thin disc 
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Figure 2. The locations of the particles with errors in vertical 
component of the gravitational acceleration colour coded as in 
Figure [1] 



analytic value. Figure |7| shows that the accuracy of this ap- 
proximation breaks down significantly in the thick disc limit. 
Note that this disc has been chosen to test the thick disc 
limit and does not represent a realistic physical situation. 
Figures [6] and [S] show a strong trend towards greater inac- 
curacy where the disc is thickest, since it is at large z that 
particles 'see' the gravitational infiuence of radial gradients 
in the disc. The accuracy also decreases at the disc's edges 
where the finite nature of disc has the largest effect. 

Figures [5] and [lO] test the accuracy of the total sur- 
face density map for the thin and thick discs, i.e. we here 
compare the result of interpolating the column density map 
obtained by sampled measurements of the projected parti- 
cle distribution with the analytic value. It is clear that the 
accuracy is relatively insensitive to disc thickness and that 
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Parameter 

H/R 

Ro-Rj 
Ro 

N 



Thin Disc 



Thick Disc Smah Disc High resolution 



0.95 0.95 
10^ 10^ 



0.02by^R/Ri 
0.33 
10^ 



0.025y^R/Ri 
0.95 
5*10^ 



Table 1. The relevant properties for the gravitational force and column density, for different analytic disc geometries. H/R is the aspect 
ratio of the disc, given as a function of cylindrical radius, Ro and Ri are the disc's outer and inner radii and N is the number of particles 
in the simulation. All discs are run with q = / Mdisc = 0.1. 
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Figure 3. The log2 ratio of the vertical component of the oct- 
tree and continuous gravitational accelerations, as a function of 
cylindrical radius for a subset of 5000 particles in the 'thick disc' 
calculation {H/R = 0.76 R/Ri). The boxplot represents the 
distribution of errors on this plot (see Figure^ for explanation). 



Spatial distribution of errors due„ 
to quantization (Thick disc) 




Figure 4. The locations of the particles with errors in vertical 
component of the gravitational acceleration colour coded as in 
Figure [3] 



most points are accurate to within a few 10s of %, with a 
median accuracy of ~ 10% for both discs. 

Overall, we find that the accuracy of our estimation 
of the column density to the surface is controlled by errors 
in determining the column density to the mid-plane, rather 
than by errors in the total surface density. The main source 
of the error in determining the column density to the mid- 



Figure 5. The log2 ratio of the vertical component of the grav- 
itational acceleration calculated using the infinite slab approxi- 
mation to the "continuous" value of the same quantity is plotted 
as a function of cylindrical radius for a subset of 5000 particles in 
the 'thin disc' calculation {H/R = 0m5y^R/Ri). The boxplot 
represents the distribution of errors on this plot (see Figure [l] for 
explanation). 



Spatial distribution of errors due to s^l^b 
approximation (Thin disc) " 
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plane depends on disc thickness (and resolution, see 3.2) 



Figure 6. The locations of the particles with errors in the vertical 
component of the gravitational acceleration due to the infinite 
slab approximation colour coded as in Figure [T] 



thin discs are dominated by discreteness effects whereas in 
thick discs the accuracy of the infinite slab approximation 
is the limiting factor. 



© 0000 RAS, MNRAS 000, 000-000 



6 Young et al. 




R/Ri R/R, 



Figure 7. The log2 ratio of the vertical component of the grav- 
itational acceleration calculated using the infinite slab approxi- 
mation to the "continuous" value of the same quantity is plotted 
as a function of cylindrical radius for a subset of 5000 particles 
in the 'thick disc' calculation (H/R = 0.76 R/Ri). The boxplot 
represents the distribution of errors on this plot (see Figure^ for 
explanation). 



Spatial distribution of errors due to slab 
approximation (Thick disc) I ^'^ 




R/R, 

Figure 8. The locations of the particles with errors in vertical 
component of the gravitational acceleration colour coded as in 
Figure [7| 



Figure 9. The log2 ratio of the calculated total surface density 
map to the known analytic values is shown as a function of radius 
for 5000 particles from the 'thin disc' simulation. The boxplot 
represents the distribution of errors on this plot (see Figure [l] for 
explanation). 




Figure 10. The log2 ratio of the calculated total surface density 
map to the known analytic values is shown as a function of radius 
for 5000 particles from the 'thick disc' simulation. See Figure 
The boxplot represents the distribution of errors on this plot (see 
Figure^ for explanation). 



3.2 Effect of resolution 

On theoretical grounds, we expect that the error due to 
quantisation will roughly depend on the "extra" accelera- 
tion imparted on a particle by its neighbours. Given that 
all particles are identical, the mass of each particle is given 
by m = Mdisc/N, where N is the number of particles in 
the simulation. Furthermore, the mean particle separation, 
A, will scale roughly as A oc N~'^^^ in a three-dimensional 
environment. Therefore, we expect the inaccuracy due to 
quantization to scale as: 

Gm 1 

This effect will be mitigated by the inclusion of gravi- 
tational softening, which for a typical SPH simulation will 



modify gravity so the effect of the quantized particles be- 
comes: 

^"^ (12) 



where s is typically set to the SPH smoothing length h and 
h is given by. 



u3 3 

pn — V m 



(13) 



where zy is a numerical parameter usually set to 1.2. Com- 
bining equations [12] and [13] we see that the inclusion of grav- 
itational softening does not change the N~'^^^ dependence 
of the quantization error, it merely decreases it by a factor 
of 1/(1 + z/)^ ^ 0.25, compared to no gravitational softening. 

In order to validate this predicted dependence on resolu- 
tion, we repeat the thin disc simulation with 5 times as many 



© 0000 RAS, MNRAS 000, 000-000 



Gravity estimate of Column Density 7 



Effect of resolution on errors 
due to quantization 
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Figure 11. The log2 ratio of the vertical component of the oct- 
tree to the continuous gravitational acceleration for a subset of 
5000 particles as a function of cylindrical radius in the 'thin disc' 
calculation ( H/R = 0.025 R/Ri) at two different resolutions 
(10^ and 5 * 10^ particles). The left and right boxplots show the 
low and high resolution data respectively. The boxplots represents 
the distribution of errors on this plot (see Figure [l] for explana- 
tion). 



n 

We 



particles (5 * 10^). The two runs are compared in Figure 
which shows that inaccuracy decreases as N increases, 
calculate the median relative error for both the high and low 
resolution simulations and find they have a ratio of ^ 1.6. 
Comparing this with equation [TT] which predicts a decrease 
in error ^ 5^^^ = 1.7, we see that our simulations are consis- 
tent with our theoretical predictions. All other comparisons 
have not been shown as they are essentially unchanged with 
resolution. 

We see that for both our N = 10^ and AT = 5 * 10^ 
runs, the error in column density estimate to the mid-plane 
is dominated for most particles by the inaccuracy of the infi- 
nite slab approximation (compare boxplots in FigurepT]with 
that in Figure [sj rather than by (resolution dependent) er- 
rors in computing Qz. Considering the region at R/Ri ^ 16 
(where H/R ^0.1 in the thin disc, a value typical of proto- 
stellar discs) we can estimate (using the above N~'^^^ scal- 
ing of discreteness errors) that N would have to be reduced 
to of order 10^ before discreteness errors dominated. Such 
a low value of N would never be employed in a hydrody- 
namic calculation in any case, since it would imply that the 
smoothing length was larger than the disc's vertical scale 



height (see Lodato & Clarke (2011) and references therein). 
We therefore conclude that our method of estimating the 
column density to the mid-plane imposes no extra resolu- 
tion requirements on the code. 



3.3 Other parameters: disc mass and radial extent 

The accuracy of our method is independent of the ratio of 
the disc mass to central object mass since the gravitational 
contribution due to the star is not used by our method. 
Changing the radial extent of the disc means that edge ef- 
fects affect the accuracy of the infinite slab approximation 
- for example, we found that a radially restricted disc (with 
fractional width of ^ 0.33 compared with ^ 0.95 in the sim- 



ulations discussed above) increased such errors by a factor 
of a few. 

Generally, the infinite slab approximation is very accu- 
rate for most particles in the simulation, with 90% of parti- 
cles being accurate to within a factor of 2 and over 50% of 
particles accurate to within a few 10s of percent. Although 
the inaccuracy is acceptable, we find that the infinite slab 
approximation slightly over-estimates the gravitational ac- 
celeration in a systematic way. Furthermore, the approxima- 
tion does worst at the edges of the disc (both radially and 
vertically) , where the deviation from an infinite slab is most 
obvious. 

The accuracy of the total surface density map produced 
by our method is high, higher even than the accuracy of the 
infinite slab approximation, which is already very good. As 
such, we do not expect the total surface density map to be 
a significant source of error in any application. 



3.4 Accuracy of the column density to the surface 
estimate. 

Having dissected the various sources in error in calculat- 
ing total surface densities and column densities to the mid- 
plane, we can understand the over-all accuracy of our de- 
termination of XI (the column density from each particle to 
the surface) which is to be used in the calculation of optical 
depths and cooling rates. Figures [12] and [13] show the rela- 
tive error in S (the column density from the particle to the 
surface), and the spatial distribution of these errors. 

It is evident that the over-all accuracy is good (most 
particles are accurate to within a factor 2) and that the 
largest errors are found at large z. This is to be expected 
since for particles close to the mid-plane, the column den- 
sity to the surface is controlled by the total surface density 
whose accuracy is of order 10%, regardless of the accuracy 
of determination of the column density to the mid-plane. At 
high the column density to the surface is obtained by the 
subtraction of two quantities of similar magnitude and the 
result is particularly sensitive to errors in determining the 
column density to the mid-plane (which, at large z, derive 
from the breakdown of the infinite slab approximation). 

In practice, however, inaccuracies at high z may be ir- 
relevant to the computation of the cooling rate according to 
equation [l] This is because once the disc enters the optically 
thin regime, the cooling rate accordingto equation [l] is in 
any case independent of optical depth. [^ 



4 APPLICATION TO REALISTIC DISCS 

Up until now we have focused on evaluating our approxi- 
mation in discs for which the column density could be an- 
alytically determined, rather than those that represent the 
outcome of realistic disc evolution. If our method is to be of 



^ This is of course not to say that equation [T] is necessarily a 
good approximation to cooling in the optically thin surface layers 
of the disc: indeed flux limited diffusion is also a poor description 
of the local cooling in such a layer and neither equation ^ nor 
flux-limited diffusion should be used in applications where the 
cooling of surface layers is of particular interest. 
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Figure 12. The relative error in E (the column density to the 
surface from a particle) is plotted for a subset of 5000 particles 
as a function of cylindrical radius in the 'thin disc' calculation 
(H/R = 0.025^^ R/Ri). The boxplot represents the distribution 
of errors on this plot (see Figure^ for explanation). 



Figure 14. The log2 ratio of (the column density between 
a particle and the mid-plane) calculated using our method and 
by direct counting for a subset of 5000 particles as a function of 
cylindrical radius for the 'marginally stable' disc (Simulation 1, 
( [Forgan et al.|20TT] >). The boxplot represents the distribution of 
errors on this plot (see Figure IT] for explanation). 



Spatial distribution of errors in E 




Figure 13. The locations of the particles with errors in E (the 
column density to the surface from a particle) colour coded as in 
Figure [121 



any use, it is important that it is accurate not only in ide- 
alised problems, but in "real world" simulations. To provide 
such a test, we evaluate our method on a disc simulation 
kindly provided by Ken Rice and detailed in Section 3 of 
(Forgan et aL|2011| (Simulation 1, Table 1). Briefly, the disc 
was evolved using SPH with radiative cooling implemented 
using the method describe in ( [Forgan et al.|200"9p . The disc 
was constructed using 5 * 10^ particles distributed between 
10 AU and 50 AU with a surface density profile E oc R~^^'^ . 
The disc to star mass ratio, was set to 0.25. The simula- 
tion was run for 27 outer rotation period^ which was long 
enough for the disc to develop spiral structures and settle 
into marginal stability (Q ^ 1), where radiative cooling was 



matched by heat generated through gravitational instabili- 
ties. 

Unlike the analytic discs, we cannot evaluate the effects 
of quantization and the infinite slab approximation sepa- 
rately. Therefore, we test our method by comparing our es- 
timated column densities to the mid-plane and the surface, 
to those obtained using a counting method described in Ap- 
pendix [C] 

Figure [M] shows a comparison of the gravity based es- 
timates of the column density to the mid-plane and the re- 
sults of the counting method. Even though this comparison 
includes inaccuracies from both the infinite slab approxima- 
tion and the quantization of the disc, the overall accuracy 
remains extremely high. In fact, over 90% of particles agree 
to within a factor of two, with most particles agreeing to 
within a few 10s of percent. Furthermore, Figure [15] shows 
that those particles with high inaccuracy are located in sur- 
face layers where they are unlikely to affect calculations of 
the cooling rate. 

Figure [16] shows the accuracy of the column density to 
the surface, E, obtained by subtracting E^ from Emap. As 
expected, the use of our total surface density map to infer 
column density to the surface does not significantly increase 
the error. The spatial distribution of the errors, shown in 
Figure [TT] show that the inaccuracy in our method only 
becomes significant in the optically thin region, where it is 
unimportant for our calculation of the cooling rate. 

For comparison, the column densities to the surface es- 
timated by the Stamatellos et al ( Stamatellos et aL]|2007 ) 



The discrep- 



method are also shown in black in Figure 
ancy is less large than that found by Wilkins & Clarke (2 012| 
(who used a thinner disc), but is significant nonetheless. 
Specifically, the Stamatellos method has a larger dispersion 
than our method and has a systematic offset of around a 



^ An outer rotation period is defined as the rotation period at 
the initial outer radius of the disc, which here is 50 AU, leading 
to an outer rotation period of 354 yrs. 



^ Only the disc's potential was used in calculating the Stamatel- 
los estimate. 
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Figure 15. The locations of the particles with errors in colour 
coded as in Figure [M] 



Figure 17. The locations of the particles with errors in E colour 
coded as in Figure [Tg] 



Errors in estimation of E 
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Figure 16. The log2 ratio of E (the column density between 
a particle and the surface) calculated using our method and by 
direct counting, for a subset of 5000 particles as a function of 
cylindrical radius for the 'marginally stable' disc (Simulation 1, 
( [Forgan et aL||2011| >). The black crosses show the log2 ratio of 
estimating E using the Stamatellos et al method ^Stamatellosj 
|et al.|[2'007| and the same value calculated by direct counting. 
The boxplot represents the distribution of errors on this plot (see 
Figure [l] for explanation). 



factor of 4. As can be seen in equation [6] the cooling rate 
depends on E quadratically in the optically thick limit, so 
this factor 4 discrepancy gives an order of magnitude under- 
estimate of the cooling rate. 

For simulations involving spiral structures, such as the 
one considered here, it is important to be able to resolve 
the differences between the cooling rates in the spiral arms 
and outside of them. In Figure [18] we plot our estimated 
total surface density map next to an "exact map" and show 
that the essential features of the disc are still clearly visible. 
This implies that our method would be able to correctly 
distinguish between on arm and off arm cooling. We find 
that (5E/E < 0.05 for < 20 and then rises steadily to a 
mean value of ^ 0.2. Comparing this with Figure |18| we 
conclude that we can easily resolve density perturbation for 



which ^E/E > 0.05. Resolving smaller density perturbations 
could be achieved by increasing from 10% the number of 
points at which the surface density map is calculateqj 



4.1 Fragmented discs 

While the above section shows that our method works well 
for self-gravitating discs in marginal stability (where heating 
due to gravitational turbulence is balanced by cooling), it 
does not test its ability to recover the cooling of fragments 
if they form. To test this, we consider a fragmented disc, 
which is originally 50 AU in size, with a disc to star mass 
ratio of .1, which has been evolved for 5 outer rotational 



periods (data provided by Peter Cossins). Figure 19 shows 
a total surface density map calculated using our method, 
which clearly shows the presence of several fragments. 

To test our methods ability to recover the column den- 
sity within a fragment, we focus the rest of our analysis on 
the .25 AU around the fragment at (-40,15) (see Figure 19). 
We calculate the column density to the surface (E) for each 
particle in this region and compare it to the column density 
obtained using the same counting method as used in section 
4] above. The results of this comparison are shown in Figure 
20] As with Figure [16] we also include the column density 
estimated by the Stamatellos method for comparison (the 
black dots in the figure). Figure 20 shows that our method 
does not perform as well as the Stamatellos method in re- 
producing the column density within fragments. Although 
our method is still accurate to within a factor of two for 
the majority of the particles, there is a systematic trend for 
our method to over-estimate the column density, while the 
Stamatellos method shows no systematically bias here, un- 
like the rest of the disc. Figure [21] shows that as before, the 
lowest accuracy points are still located at high z, where the 
column density is less important to estimating the cooling 
rate. 



^ Note the choice of 10% of particles is somewhat arbitrary. Spiral 
structure can still be resolved even if the column density is calcu- 
lated at only 1% of particles (data not shown). This is discussed 
in more detail in Appendix [B] 
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Figure 18. The two heat maps show the total surface density 
inferred by the method used to estimate the column density (top) 
and an "exact" counting method (bottom). The colour bars show 
the mapping between total surface density and colour used in this 
plot (log scale). 
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Figure 19. A heat map showing the total surface density esti- 
mate for the fragmented disc described in section [XT] The colour 
bar shows the mapping between total surface density and colour 
on a log scale. 
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Figure 20. The log2 ratio of E (the column density between 
a particle and the surface) calculated using our method and by 
direct counting, for all particles within .25 AU of the fragment 
located at (-40,15) in Figure [To] as a function of cylindrical ra- 
dius from the centre of the fragment. The black crosses show the 
log2 ratio of estimating E using the Stamatellos et al method 
( [Stamate llos et al.|2007|> and the same value calculated by direct 
counting. The boxplot represents the distribution of errors on this 
plot (see Figure [l] for explanation). 
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Figure 21. The locations of the particles with errors in E colour 
coded as in Figure [20| as a function of cylindrical radius Rfrag 
and z from the centre of the fragment. 



It is perhaps unsurprising that our method performs 
less well within a fragment, as fragments have a roughly 
spherical geometry where the infinite slab approximation is 
likely to break down. The same reasoning explains why the 
Stamatellos method does so well here, as it is known to 
perform well for problems with a spherical geometry. 

The calculation of our total surface density map pro- 
vides a mechanism for combining the advantages of our 
method with those of the Stamatellos method in fragmented 
discs. That is, the Stamatellos method can be used to esti- 
mate E within the fragments and our method can be used 
everywhere else, where a fragment can be defined to be a 
region where the total surface density exceeds the average 
value by two orders of magnitude. Indeed, this is the ex- 
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act definition used to identify fragments in many studies of 
fragmenting discs (e.g. Forgan et al. (2011)). Such a scheme 
would combine the advantages of the two methods and pro- 
vide an estimate of the column density without systematic 
bias for the entirety of the fragmented disc and may be op- 
timal when the cooling within fragments is an important 
consideration. 
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5 CONCLUSION 

In this paper we have presented a new method for efficiently 
and accurately estimating the cooling in disc geometries. 
To achieve this, our method estimates the column density 
between each particle and the surface of the disc, via esti- 
mates of the column density to the mid-plane and the total 
surface density. We verify the accuracy of our column den- 
sity estimation (and hence our cooling rates) by comparison 
with discs for which the analytic form of the column den- 
sity can be calculated. We test our method on a realistic 
proto-planetary disc simulation, that has been evolved for 
long enough to reach marginal stability and develop a typ- 
ical spiral structure. Finally, we test our method on a frag- 
mented disc and find that our method does not perform as 
well as the Stamatellos method within the fragments, due to 
the locally spherical geometry. We suggest that this short- 
coming can be resolved by using the Stamatellos method 
only within regions of extremely high density (i.e. the frag- 
ments). We find throughout our tests that the accuracy of 
our method remains high (i.e. typical errors of order a few 
tens of per cent) and conclude that it is ideally suited for 
use in problems that depend on an accurate estimate of the 
cooling rate in disc geometries. 



6 MATERIALS & METHODS 

In the interests of reproducibility and transparency, all code 
and data used in performing this work have been made freely 
available online. 

The generation of initial conditions and the code 
used to perform the analyses described in this paper can 
be found at https : / /bitbucket . org/ constant Amateur/ 
jdisccolumndensity See the readme file in this repos- 
itory for further details. The gravitational acceleration 
for the analytic discs was calculated using a modified 
version of GADGET-2.0 (Springel^ 2005) , which can be 
obtained from https : / /bitbucket . org/ constant Amateur/ 
gadgetoutputgravaccel Converting of initial conditions 
to/from ascii files was done using the code from https: 
[//bitbucket . org/ constantAmateur/ easyic 

The data from the spiral simulation used in sec- 
tion |4] were provided by Ken Rice (Forgan et al. 2011). 
The fragmented disc used in section~ |4.1| were provided 
by Peter Cossins and Giuseppe Lodato (unpublished). 
Both are made available at https : //bitbucket . org/ 
const ant Amateur /disccolumndensity with the original au- 
thors' permission. 



APPENDIX A: THE ANALYTIC COLUMN 
DENSITY 

Consider a disc with density profile {p{R, z)) given by: 



27T (^Ro' 



+2 



(Al) 



This corresponds to a disc of mass M^, inner and outer 
radii Ri and Ro and surface density profile S oc R^p; the 
Gaussian distribution with respect to z corresponds to a 
situation of hydrostatic equilibrium in the case that the disc 
is locally vertically isothermal, in which case the scale height 
H is given by: 



H 



liniH 
R^ 



(A2) 
(A3) 

(A4) 



where 7 = 5/3 is the adiabatic index, /x = 2.3 is the molecu- 
lar weight, rriH is the mass of Hydrogen and kh is the Boltz- 
mann constant. Finally, the temperature profile is given by 



Putting this all together gives: 



H = 



GM.firriHRT 



7R~ 



(A5) 



(A6) 



The vertical component of gravitational acceleration at 
R,z is given by the following integral: 

GR'p{R', z'){z' - z)dR'de'dz 
(i?2 + R'^ - 2RR' cos{0 - 0') + {z' - zYf'^ 

Since this integral cannot be solved analytically, we 
evaluated it numerically using the "NIntegrate" function in 
mathematica using the "Adaptive Monte Carlo" method. 
The column density from each particle at R^z to the sur- 
face and to the mid-plane are obtained from integration of 
equation I Al I i.e. 



(A7) 



E' 



Md(ap + 2) 



4^ (i?^+'-i?r^+') 

Mi{ap + 2) 
4^ (rZ"^^ - fir-+') 



R°'''erf 



R^ferfc 




(A8) 



(A9) 
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Md{ap + 2) 



4^ 



(AlO) 



where erf is the error function and erfc is the comple- 
mentary error function. 



APPENDIX B: CALCULATING THE TOTAL 
SURFACE DENSITY MAP 

For our method to be useful, we have to be able to convert 
column densities to the mid-plane to column densities to 
the surface without significant loss of accuracy. The most 
obvious way to do this is to try and use the fact that we know 
the column density to the mid-plane for all particles and use 
the particles at the top of the disc to approximate the surface 
density. However, as is shown in the main text, the accuracy 
of the column density to the mid-plane estimates is lowest at 
the top of the disc. As each estimate of the surface density 
will potentially effect several particles below it, even small 
inaccuracies will be compounded. 

In essence, what is required is a method for estimating 
the two-dimensional surface density from the particle po- 
sitions and masses projected onto the mid-plane. There are 
many such methods available, each with advantages and dis- 



advantages (see section 2.3 and Ferdosi et al. (2011|). Our 
method is not tied to any one technique for estimating the 
total surface density and so the advantages and disadvan- 
tages of the different techniques should be weighed against 
the scientific application of interest. In this paper, we adopt 
the following algorithm to calculate a total surface density 
map. 

(i) Randomly select 10% of the particles, call these the 
tracer particles. 

(ii) Project all particles on the xy plane and for each 
tracer particle, find its 60 closest projected neighbours 

(iii) The column density at each tracer particle is then 
given by (mass of 60 nearest particles) / TvR^ax where Rmax 
is the distance between the tracer particle and its 60th fur- 
thest neighbour 

Because the tracer particles are chosen at random, the 
resolution of our surface map automatically adjusts to the 
density profile of our disc. Furthermore, because we are cal- 
culating the surface density directly, the inaccuracies of the 
gravity based estimates at large z are not an issu^ This 
method has a computational complexity that is roughly 
OiNlogN). 

An alternative of only slightly lower accuracy (data not 
shown), but significantly improved computational efficiency 
0{N), is to construct a two dimensional grid and evaluate 
the column density by adding up the masses in each grid 
cell and dividing by the grid area. The grid can be spaced 
so that roughly the same number of particles is located in 



^ The accuracy of the estimate of surface density at (x,y) could 
be improved by using the 2D version of an SPH smoothing here. 
However, this map will be used by many points in the cylinder 
to represent the column density to the surface, so it is preferable 
that our estimate be a little more "washed out" to increase the 
accuracy for those particles off the (x,y) axis that also use this 
grid point as an estimator of E. 



each radial annulu^ allowing the grid to adjust to the disc's 
density profile. The surface density for each particle is then 
given by the cell within which it resides. 

The exact computational cost of both the "random sam- 
pling" and "grid" methods suggested here will in general 
depend upon the code which is used for the rest of the simu- 
lation (i.e. the hydrodynamics and gravity). This is because 
different codes have different costs associated with cross pro- 
cess communication and different logical times where all par- 
ticles in the simulation are easily accessible. 



APPENDIX C: CALCULATING AN "EXACT" 
COLUMN DENSITY 

For the realistic simulation considered in section [4] the col- 
umn density at each point to the surface or the mid-plane 
cannot be determined analytically. As such, we require some 
method to calculate the column density that we can use as 
our gold standard , to test the accuracy of our method. To 
do this we use the same basic idea as employed in Appendix 
|B]to build the total surface density map. In detail, we do 
the following for each particle. 

(i) Remove all particles from the simulation that are be- 
low the particle (for column density to surface) or not be- 
tween the particle and the mid-plane (for column density to 
the mid-plane). 

(ii) Project all remaining particles onto the xy plane, find 
either the 60 closest particles or the number of particles 
within a circle of radius 3 AU, centred on the point for which 
we are trying to calculate the column density. 

(iii) The column density is then given by the sum of 
the masses of our neighbouring particles divided by TrR^ax 
which Rmax is either the distance to the 60th particle (in 
the xy plane) or 3 AU if there are fewer than 60 neighbours 
within 3 AU. 
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